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Abstract.  We  develop  a  second-order  rotated  grid  method  for  the  approximation  of  time  de¬ 
pendent  solutions  of  conservation  laws  in  complex  geometry  using  an  underlying  Cartesian  grid. 
Stability  for  time  steps  adequate  for  the  regular  part  of  the  grid  is  obtained  by  increasing  the  domain 
of  dependence  of  the  numerical  method  near  the  embedded  boundary  by  constructing  /i-boxes  at 
grid  cell  interfaces.  We  describe  a  construction  of  h-boxes  that  not  only  guarantees  stability  but 
also  leads  to  an  accurate  and  conservative  approximation  of  boundary  cells  that  may  be  orders  of 
magnitude  smaller  than  regular  grid  cells.  Of  independent  interest  is  the  rotated  difference  scheme 
itself,  on  which  the  embedded  boundary  method  is  based. 
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1.  Introduction.  We  present  a  high-resolution  Cartesian  grid  embedded  bound¬ 
ary  method  for  the  approximation  of  hyperbolic  conservation  laws  in  complex  geome¬ 
tries.  A  Cartesian  grid  approach  is  attractive,  since  away  from  the  boundary  it  allows 
the  use  of  standard  high-resolution  shock  capturing  methods  that  are  more  difficult 
to  develop  on  unstructured  (body  fitted)  grids.  Embedded  boundary  grids  allow  a 
more  automated  grid  generation  procedure  around  complex  objects,  which  is  impor¬ 
tant  especially  for  three-dimensional  problems.  The  numerical  challenge  associated 
with  a  Cartesian  grid  embedded  boundary  approach  is  the  so-called  small  cell  prob¬ 
lem.  Near  the  embedded  boundary  the  grid  cells  may  be  orders  of  magnitude  smaller 
than  regular  Cartesian  grid  cells.  Since  standard  explicit  finite  volume  methods  take 
the  time  step  proportional  to  the  size  of  a  grid  cell,  this  would  typically  require  small 
time  steps  near  an  embedded  boundary.  We  will  describe  a  numerical  method  that 
overcomes  the  time  step  restriction,  while  seeking  an  accurate  approximation  near  the 
embedded  boundary  as  well  as  in  the  whole  domain. 

During  the  last  decade,  many  different  Cartesian  grid  embedded  boundary  meth¬ 
ods  have  been  developed.  Several  authors  use  a  cell  merging  technique,  where  small 
irregular  cut  cells  are  merged  together  with  a  neighboring  regular  grid  cell,  see  for 
instance  Coirier  and  Powell  [9]  or  Quirk  [33].  An  interesting  variant  of  cell  merging 
was  suggested  by  Falcovitz  et  al.  [14],  who  combined  cell  merging  with  dimension 
splitting.  Chern  and  Colella  [8]  developed  a  flux  redistribution  procedure  for  front 
tracking  methods.  This  approach  was  extended  by  Pember  et  al.  [32]  and  Modiano  and 
Colella  [29]  to  overcome  the  small  cell  problem  for  an  embedded  boundary  method. 
In  [16],  Forrer  and  Jeltsch  suggested  an  embedded  boundary  method  where  all  grid 
cells  are  handled  as  full  cells.  Grid  cells  near  the  boundary  are  updated  in  each  time 
step  to  simulate  a  reflecting  boundary.  A  similar  approach  was  also  suggested  by  Xu 
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et  al.  [38].  Although  the  small  cell  problem  is  overcome  by  all  of  these  approaches, 
the  accuracy  near  the  boundary  is  reduced.  Some  of  the  available  methods  (i.e.  [16], 
[38])  are  not  conservative,  although  this  error  may  be  relatively  small. 

Alternatively  one  can  use  an  implicit  scheme  at  the  cut  cells,  but  unless  this  is 
the  method  of  choice  for  the  rest  of  the  domain  [31],  this  is  an  expensive  alternative. 

An  accurate  approximation  near  the  boundary  may  in  particular  be  important 
for  the  simulation  of  flow  problems  with  moving  objects.  Cartesian  grid  methods  for 
conservation  laws  in  the  presence  of  moving  objects  have  for  instance  been  studied  in 
the  literature  by  Falcovitz  et  al.  [14]  (using  cell  merging)  and  Forrer  and  Berger  [15] 
(extending  the  approach  presented  in  [16]).  For  the  test  case  studied  in  those  papers 
the  two  different  methods  produced  a  different  solution  structure.  Since  the  motion 
of  the  embedded  object  depends  on  the  numerical  solution  near  the  boundary,  this 
may  indicate  a  strong  need  for  the  development  of  more  accurate  embedded  boundary 
methods. 

Our  approach  to  overcome  the  small  cell  problem  is  based  on  the  so-called  h- 
box  method  suggested  by  Berger  and  LeVeque  [4],  [5].  Here  we  will  extend  this 
method  to  a  high-resolution  scheme.  The  basic  idea  behind  the  h- box  method  is  to 
approximate  numerical  fluxes  at  the  interface  of  a  small  cell  based  on  initial  values 
specified  over  regions  of  length  h ,  where  h  is  the  length  of  a  regular  grid  cell.  In 
order  to  obtain  the  required  stability  property,  the  flux  differences  at  cut  cells  should 
be  of  the  size  of  the  small  irregular  cell.  Here  we  develop  a  high-resolution  rotated 
grid  method  that  has  this  cancellation  property  of  numerical  fluxes.  We  first  study 
our  high-resolution  rotated  grid  h- box  method  for  regular  Cartesian  grids.  Beside  the 
development  of  an  accurate  embedded  boundary  method,  the  rotated  grid  method 
that  is  described  in  this  paper  may  be  of  interest  also  for  other  applications.  A 
crucial  step  in  the  development  of  this  method  is  the  definition  of  the  left  and  right 
states  that  are  used  in  the  flux  calculation  of  the  rotated  grid  approach.  We  will  show 
that  for  the  approximation  of  the  two-dimensional  advection  equation,  our  rotated 
grid  h- box  method  is  similar  to  an  upwind  method  described  by  Roe  and  Sidilkover, 
see  [34].  This  implies  an  interesting  optimality  property.  For  the  flux  calculation 
a  predictor-corrector  approach,  based  on  Colella’s  multidimensional  high-resolution 
method  described  in  [10],  is  used.  We  will  show  that  in  the  context  of  a  rotated  grid 
method  the  predictor-corrector  approach  has  advantages  over  other  existing  methods. 

Several  different  rotated  grid  methods  for  the  approximation  of  multi-dimensional 
systems  of  conservation  laws  can  be  found  in  the  literature,  see  for  instance  [13],  [18], 
[27],  [34],  The  main  motivation  for  developing  these  methods  was  to  reduce  grid  ef¬ 
fects  by  calculating  fluxes  in  a  direction  that  depends  on  the  solution  structure  instead 
of  the  grid.  In  [27],  it  was  shown  that  rotated  grid  methods  can  lead  to  a  sharper 
approximation  of  shock  waves  that  are  not  aligned  with  the  grid.  In  [26],  LeVeque 
and  Walder  explored  the  use  of  first  order  accurate  rotated  grid  h- box  methods  for 
the  approximation  of  astrophysical  applications.  Nevertheless,  what  is  really  needed 
for  the  approximation  of  conservation  laws  are  high-resolution  methods  that  are  sec¬ 
ond  order  accurate  in  smooth  regions  and  lead  to  a  sharp  approximation  of  shock 
waves  while  spurious  oscillations  are  suppressed.  Compared  to  the  rotated  grid  h- box 
method  used  in  [26]  as  well  as  [4],  [5],  our  new  method  uses  a  different  h- box  length, 
another  definition  of  h- box  values  as  well  as  a  different  flux  calculation.  All  of  those 
changes  were  necessary  in  order  to  obtain  a  method  that  is  second  order  accurate  for 
smooth  solutions  as  well  as  stable  for  time  steps  with  a  Courant  number  of  up  to  one. 

In  [2],  we  studied  the  construction  of  high-resolution  h- box  methods  for  the  ap- 
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proximation  of  one-dimensional  systems  of  conservation  laws  on  irregular  grids,  see 
also  [6]  for  earlier  related  work.  This  provides  basic  insight  for  the  construction  of  an 
embedded  boundary  method  in  multi-dimensions.  For  the  simpler  one-dimensional 
situation  we  showed  that  our  h- box  method  leads  to  a  stable  and  second  order  accu¬ 
rate  approximation  of  the  advection  equation  on  arbitrary  irregular  grids.  Numerical 
results  confirmed  the  same  accuracy  and  stability  properties  for  the  approximation 
of  the  Euler  equations.  We  are  not  familiar  with  other  finite  volume  methods  for 
one  dimensional  conservation  laws  on  irregular  grids  that  have  the  same  accuracy  and 
stability  property.  The  large  time  step  Godunov  method  of  LeVeque  described  in  [20], 
[21],  [22]  is  related  to  the  h- box  method.  This  scheme  allows  larger  time  steps  in  the 
approximation  of  nonlinear  systems  of  conservation  laws  by  increasing  the  domain  of 
influence  of  the  numerical  scheme.  This  is  done  in  a  wave  propagation  approach,  in 
which  waves  are  allowed  to  move  through  more  than  one  mesh  cell.  The  interaction  of 
waves  is  approximated  by  linear  superposition.  At  a  reflecting  boundary  this  method 
becomes  more  difficult  than  an  h- box  method  especially  in  higher  dimensions,  since 
the  reflection  of  waves  at  the  boundary  has  to  be  considered  for  waves  generated  by 
Riemann  problems  away  from  the  boundary,  see  [3],  In  [30],  Lemma  3.5,  Morton 
showed  that  high-resolution  versions  of  such  a  large  time  step  method  lead  to  a  sec¬ 
ond  order  accurate  approximation  of  the  one-dimensional  advection  equation  on  a 
non-uniform  grid  only  if  the  grid  varies  smoothly.  Our  second  order  h- box  method 
does  not  require  this  smoothness  assumption. 

The  rest  of  the  paper  is  organized  as  follows.  In  section  2,  we  develop  the  high- 
resolution  rotated  grid  h- box  method  for  the  approximation  of  conservation  laws  on 
regular  two-dimensional  Cartesian  grids.  It  is  a  multidimensional  method,  in  that  it 
does  not  use  dimension  splitting.  Our  new  embedded  boundary  rotated  grid  h- box 
method  will  be  described  in  section  3.  Throughout  this  paper  we  use  the  advection 
equation  (with  constant  as  well  as  varying  velocity  field)  and  the  Euler  equations  of 
gas  dynamics  to  illustrate  the  numerical  method.  Numerous  simulations  demonstrate 
the  performance  of  our  approach. 

2.  A  high-resolution  rotated  grid  method  for  regular  Cartesian  grids. 

We  consider  the  approximation  of  a  two  dimensional  system  of  conservation  laws  on 
a  regular  Cartesian  grid.  Under  appropriate  smoothness  assumptions  the  equations 
can  be  written  in  the  differential  form 

^ f(q{x,y,t ))  +  g(q{x,y,t ))  =  0,  (2.1) 

where  q(x,  y,  t )  is  a  vector  of  conserved  quantities  and  f(q(x ,  y,  t))  and  g(q(x,  y ,  t))  are 
vector  valued  flux  functions.  For  the  numerical  approximation,  we  use  a  conservative 
finite  volume  method  that  takes  the  general  form 


The  numerical  fluxes  F  and  G  at  the  cell  interfaces  are  calculated  via  a  rotated  grid 
approach  with  coordinate  directions  £  and  rj  and  angle  of  rotation  6.  In  the  new 
coordinate  system  the  hyperbolic  equation  (2.1)  takes  the  form 
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with 

/(<?)  =  /(<?)  cos (0)  +  g(q)  sin(0) ,  (2.4) 

§(q)  =  ~f(q)  sin(0)  +  5(g)  cos(0).  (2.5) 

At  each  cell  interface  of  the  Cartesian  grid  we  solve  Equations  (2.4)  and  (2.5)  numer¬ 
ically.  This  gives  the  numerical  fluxes  F  and  G.  The  numerical  fluxes  F  and  G  that 
are  used  in  the  finite  volume  scheme  (2.2)  can  then  be  obtained  via  the  relations 

F  —  Fcos(9)  -  Gsin(0)  (2.6) 

G  —  F  sin(0)  +  G  cos(0) .  (2.7) 

The  numerical  fluxes  F  and  G  are  calculated  by  solving  Riemann  problems  for  the 
one-dimensional  equations 


|?  +  |/(?)  =  0 

(2.8) 

d  d  ,  , 

(2.9) 

^+^5(9)=°. 

The  crucial  step  in  order  to  obtain  an  accurate  high-resolution  rotated  grid 
method  is  the  definition  of  left  and  right  states  of  the  conserved  quantities  that  are 
used  as  initial  values  for  the  solution  of  Riemann  problems  in  the  £  and  ^-direction 
respectively.  We  will  derive  a  new  rotated  grid  h- box  method  that  is  motivated  by  a 
multidimensional  upwind  method  of  Roe  and  Sidilkover,  see  [34],  For  simplicity  we 
first  restrict  our  consideration  to  the  approximation  of  the  two-dimensional  advection 
equation. 

2.1.  A  rotated  grid  method  for  advective  transport.  Consider  the  two- 
dimensional  advection  equation  with  constant  advection  speed,  i.e. 

Q  0  0 

—q(x,y,t)  +  u—q(x,y,t)  +  v—q(x,y,t)  =  0.  (2.10) 

In  order  to  develop  a  high-resolution  rotated  grid  method  we  want  to  use  appro¬ 
priate  linear  interpolation  formulas  to  define  the  conservative  quantities  used  in  the 
flux  calculations.  Note  that  this  was  also  the  crucial  step  in  developing  high-resolution 
irregular  grid  h- box  methods  in  the  ID  case,  see  [2].  We  will  formulate  the  rotated 
grid  method  as  an  h- box  method.  The  use  of  h- boxes  allows  simple  geometric  inter¬ 
pretations  which  will  later  be  useful  for  the  construction  of  the  embedded  boundary 
method.  To  make  comparisons  with  other  existing  schemes  easier,  we  first  restrict 
our  considerations  to  the  case  0.5  <  u/v  <  2.  We  will  see  below  that  our  rotated  grid 
h- box  method  generalizes  in  a  straightforward  way  to  arbitrary  values  of  the  advection 
speed.  Our  approach  for  finding  the  h- box  values  is  motivated  by  multidimensional 
upwind  methods  described  in  Roe  and  Sidilkover  [34],  where  the  approximation  of 
steady  state  solutions  for  the  advection  equation  was  considered.  See  also  [36]  and 
[37]  for  related  work. 

A  geometric  interpretation  of  a  two-dimensional  second  order  accurate  upwind 
method  for  steady  state  problems  is  shown  in  Figure  2.1a.  At  each  cell  interface  the 
conservative  quantity  that  is  used  in  the  upwind  scheme  is  defined  by  tracing  back  the 
characteristic  at  the  midpoint  of  the  cell  interface  (for  instance  A).  Linear  interpola¬ 
tion  onto  the  Cartesian  grid  using  C  and  D  is  used  to  find  the  values  that  are  assigned 
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Fig.  2.1.  Schematic  description  of  the  advection  algorithm  suggested  by  Roe  and  Sidilkover 
[3^].  (cl)  A  second  order  residual  results  from  interpolating  to  the  values  at  the  open  squares,  (b) 
The  optimum  first  order  method  (N-scheme)  results  from  moving  each  interpolation  point  to  the 
midpoint  of  its  interval,  if  initially  it  is  more  than  half  a  length  away  from  the  center  of  the  grid 
cell  interface  considered. 


to  the  characteristic.  In  [36,  section  4.2.1]  it  was  shown  that  the  resulting  upwind 
method  is  second  order  accurate  for  smooth  steady-state  solutions  of  the  advection 
equation.  However  this  scheme  is  not  monotone  and  therefore  not  suitable  for  the 
approximation  of  discontinuous  solutions.  In  [34],  Roe  and  Sidilkover  derived  an  op¬ 
timal  positive  method,  i.e.  a  first  order  accurate  method  with  the  smallest  truncation 
error  among  all  monotone  schemes  that  use  the  same  narrow  stencil.  A  geometric 
interpretation  of  this  method  is  illustrated  in  Figure  2.1b.  The  upwind  points  for  the 
flux  calculations  Gij+i  and  G%]  1  are  obtained  by  characteristic  interpolation  while 
the  upwind  points  used  for  the  fluxes  Fi  1  j  and  Fi+l  j  are  obtained  by  arithmetic  av¬ 
eraging,  i.e.  the  interpolation  points  for  the  calculation  of  the  fluxes  in  the  ^-direction 
are  moved  upward.  Assuming  that  the  advection  speeds  satisfy  0.5  <  u/v  <  2,  fluxes 
are  found  by  characteristic  interpolation,  if  the  interpolation  point  is  at  most  a  half 
grid  cell  length  apart  from  the  center  of  the  cell  interface  (i.e.  the  G  fluxes  in  Figure 
2.1a),  otherwise  the  interpolation  point  is  moved  and  the  corresponding  fluxes  are 
determined  by  averaging.  This  method  is  often  called  the  N-scheme  for  its  narrow 
stencil.  Here  we  will  discuss  an  interpretation  of  the  N-scheme  as  a  rotated  grid  h- box 
method  that  will  lead  to  a  natural  extension  to  a  high-resolution  rotated  grid  method 
for  time  dependent  problems. 

Figure  2.2  shows  an  interpretation  of  the  optimal  positive  scheme  depicted  in 
Figure  2.1b  as  an  h- box  method.  At  each  cell  interface  we  construct  h- boxes  in  the 
characteristic  direction,  i.e.  in  direction  (it,  v)  for  the  advection  equation.  In  Figure 
2.2a  we  show  h- boxes  that  are  used  for  the  calculation  of  G  fluxes.  The  interpolated 
point  (illustrated  in  Figure  2.1b)  that  was  used  for  the  approximation  of  the  G  fluxes 
by  the  optimal  positive  scheme  can  also  be  obtained  by  calculating  the  corresponding 
h- box  values  based  on  the  area  fraction  of  the  grid  cells  that  are  overlapped  by  the 
h- box.  The  N-scheme  is  obtained  if  we  set  the  length  of  the  h- boxes  equal  to 

h  —  y/l  +  (min(|u|,  \v\)/  max(|u|,  |i;|))2Ax.  (2.11) 

Note  that  in  earlier  work  h- boxes  of  length  h  —  Ax  —  Ay  were  used,  see  [4],  [5], 
[26].  For  simplicity  we  assume  Ax  —  Ay  throughout  this  paper.  With  the  /i-box 
length  (2.11)  the  interpolation  used  for  the  calculation  of  the  fluxes  F  by  the  optimal 
positive  scheme  (compare  with  Figure  2.1b)  also  has  a  geometric  interpretation.  Again 
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Fig.  2.2.  H -boxes  in  the  £- direction  for  the  calculation  of  the  fluxes  G  and  F  respectively,  that 
correspond  to  the  optimal  positive  scheme. 


averaging  over  the  h- boxes  depicted  in  Figure  2.2b  suggests  using  the  average  of  the 
conserved  quantities  in  the  two  grid  cells  that  are  overlapped  by  the  h-box.  This 
h- box  method  is  stable  for  time  steps  that  satisfy  CFL  <  1.  To  show  stability  we 
restrict  our  considerations  to  the  case  u  <  v,  with  u,  v  >  0.  The  method  can  be 
written  in  the  form  of  a  finite  volume  method 


KiJ  -  «?-!•>)  -  £«  («U  -  0«-l) '  <212) 


Ay 


The  left  h- boxes  used  in  the  above  equation  have  the  form: 


,j-i)  > 


Q 


j  —  2  (Qi-l,j 


Q 


i-l,j-l) 


Qi+1  ,j  ~  2 

QA+h  =  YvQi-^  +  (1_  2^)  Qi'j’  QA  I  =  +  V1  “  2^)  Qi’j-1' 

(2.13) 


Using  these  expressions  we  can  rewrite  the  finite  volume  method  (by  again  using 
Ax  —  Ay)  into  the  form 


QZV^Q?Ai-v^)+Qri 


**>3 


At 


Ax 


*i,3- 


At 

^V~u)Ai 


Q1 


i- l,j— 1  * 


<2'14> 


All  coefficients  on  the  right  hand  side  of  Equation  (2.14)  are  greater  than  or  equal  to 
zero  assuming  that  v-^  <  1.  In  the  general  case,  i.e.  without  the  above  assumptions 
on  u  and  v,  the  time  step  restriction  for  stability  takes  the  form  At  max(|u|,  |u|) / Ax  < 
1. 

Other  multi-dimensional  methods  for  the  advection  equation  can  also  be  inter¬ 
preted  as  rotated  grid  h- box  methods.  If  we  reduce  the  h- box  length  to 

h  —  At u2  +  v2, 


then  h- box  values  that  are  again  calculated  by  averaging  over  piecewise  constant  cell 
averages  of  the  conserved  quantity  take  the  form 


Qi,j+,  —  Qij 


vAt 
2A y 
uAt 
2  Ax 


( Qij 


( Qi,j 


Qi,j—i)’  Qj-ij  —  Qi-i,j  2a 1d 
Qij  —  i  —  Qi,j- 1  —  2Ax 


Qi-lJ-l) 

Q i  —  l,j  —  l) ' 

(2.15) 
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The  resulting  rotated  grid  h- box  method  (2.12),  (2.15)  is  equivalent  with  the  corner 
transport  upwind  method  (CTU)  of  Colella  [10].  The  CTU  method  can  also  be 
interpreted  as  wave  propagation  algorithm,  see  [23],  [24].  This  method  is  also  first 
order  accurate  and  stable  for  time  steps  that  satisfy  CFL  <  1.  In  this  paper  we  will 
concentrate  on  h- box  methods  with  the  fixed  h- box  length  given  by  Equation  (2.11). 

Note  that  only  in  the  special  cases  u  —  v  —  1  (i.e.  6  —  tv/ 4,),  u  —  {)  or  v  — 
0  does  the  optimal  positive  scheme  of  Roe  and  Sidilkover  lead  to  a  second  order 
accurate  approximation  in  space.  For  time  dependent  problems  we  want  to  obtain 
a  high-resolution  approximation  by  solving  Riemann  problems  (defined  by  the  h- box 
values)  with  a  high-resolution  method.  Note  that  by  using  a  Lax-Wendroff-type  flux 
calculation  (i.e.  to  get  2nd  order  in  time)  combined  with  the  h- box  method  that 
corresponds  to  Roe  and  Sidilkover’s  approach  we  only  get  a  second  order  accurate 
approximation  in  the  special  cases  mentioned  above. 

There  are  high-resolution  versions  of  the  N-scheme,  see  for  instance  [37].  In  order 
to  combine  the  monotonicity  property  with  second  order  accuracy  a  nonlinear  scheme 
was  constructed  which  blends  between  the  N-scheme  and  the  second  order  accurate 
non-monotone  method.  To  retain  monotonicity  limiters  are  applied  to  the  difference 
between  the  fluxes  of  the  two  methods.  However,  those  methods  cannot  be  interpreted 
as  a  rotated  grid  h- box  method.  As  we  will  see  in  section  3,  stability  of  the  embedded 
boundary  method  relies  on  the  use  of  a  rotated  grid  method  in  which  the  upwind  point 
can  be  interpreted  as  an  approximation  of  the  conserved  quantities  at  the  center  of 
mass  of  an  h- box.  Therefore,  we  will  now  develop  a  new  second  order  accurate  rotated 
grid  h- box  method  that  has  this  property. 

In  order  to  obtain  a  second  order  accurate  approximation  with  a  rotated  grid 
h- box  method,  the  h- box  values  have  to  be  calculated  by  appropriate  linear  inter¬ 
polation.  Our  aim  is  to  construct  a  high-resolution  rotated  grid  h- box  method  that 
leads  to  second  order  accurate  results  for  smooth  solutions  and  avoids  spurious  os¬ 
cillations  for  the  approximation  of  discontinuities.  We  will  derive  such  a  scheme  by 
using  slope  limiters  in  our  interpolation  scheme.  If  all  slopes  are  set  to  zero,  our 
scheme  will  reduce  to  the  optimal  positive  scheme  of  Roe  and  Sidilkover.  Figure  2.2b 
shows  a  description  of  the  interpolation  that  is  used  to  obtain  the  h- box  values.  We 
calculate  the  conserved  quantities  assigned  to  the  left  h- box  by  first  determining  the 
conserved  quantities  at  the  points  I\  and  I 2  shown  in  the  figure.  Those  values  are 
in  the  unlimited  case  obtained  by  linear  interpolation  between  neighboring  grid  cell 
values  (e.g.  QtJ  and  Qi+ij  are  used  to  calculate  the  conserved  quantities  at  I2).  The 
conservative  quantities  assigned  to  the  h- box  are  then  obtained  as  the  arithmetic  av¬ 
erage  between  the  conserved  quantities  at  the  points  Ji  and  I2.  Note  that  the  h- box 
values  can  be  considered  as  a  linear  interpolation  at  the  center  of  mass  of  the  h- box 
(i.e.  the  little  filled  squares  marked  in  Figure  2.2  are  located  at  the  center  of  mass 
of  the  h- box).  In  the  case  shown  in  Figure  2.2a,  simple  averaging  over  the  region 
overlapped  by  the  h- box  is  equivalent  to  linear  interpolation  to  the  center  of  mass  of 
the  h- box.  The  resulting  h- box  values  can  be  used  to  construct  second  order  accurate 
methods  with  appropriate  second  order  flux  formulas.  In  this  case  no  slope  limiting 
in  the  interpolation  step  for  the  calculation  of  the  h- box  values  is  necessary. 

To  avoid  unphysical  oscillations  near  discontinuities  (even  for  the  first  order  up¬ 
date)  we  include  slope  limiters  into  the  interpolation  process  used  for  h- box  values 
of  the  form  indicated  in  Figure  2.2b.  This  is  described  in  the  next  section.  Once 
the  h- box  values  are  defined  we  can  calculate  the  fluxes  in  the  £  and  //-directions  by 
a  first  order  or  a  high-resolution  method.  This  is  analogous  to  the  one-dimensional 
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irregular  grid  h- box  method  described  in  [2],  To  obtain  flux  limiters  that  are  needed 
by  a  high-resolution  flux  calculation,  we  construct  an  additional  h- box  on  each  side. 

2.2.  Slope  limiters  for  the  interpolation  process.  In  this  section  we  will 
restrict  our  considerations  to  the  approximation  of  the  advection  equation  with  posi¬ 
tive  advection  speed  u ,  v  and  assume  v  >  u.  We  choose  our  ^-direction  of  the  rotated 
grid  method  to  be  aligned  with  the  constant  velocity  field.  For  the  advection  equation 
this  implies  that  we  do  not  have  to  calculate  flux  components  in  the  ^-direction.  Our 
interpolation  method  for  the  calculation  of  the  h- box  values  requires  the  use  of  slope 
limiters  to  avoid  unphysical  oscillations  near  discontinuities.  Such  a  slope  limiter  is 
already  necessary  even  if  the  fluxes  are  calculated  with  the  first  order  upwind  method, 
and  we  therefore  restrict  our  considerations  in  this  subsection  to  this  case. 

Our  numerical  method  can  be  written  in  the  conservative  form  (2.2),  with 

Fk+i,j  =uQ%+ij,  ke{i-l,ij 
Gi,k+i  =  vQt,k+p  ke{j-l,j} 

The  left  h- box  values  are 

Qk+bi  ~  2  +  Qkd)  +  2  +  Sk+bi)  dl'  dl  ~  2  _  2v'  k  ^  ~ 

Qi,k+i  —  (1  “  d2)Qi,k  +  d2Qi-i,k,  d2  —  — ,  k  e  {j  -  1  ,j}. 

(2.16) 

In  the  above  formulas  di  is  the  distance  of  the  interpolated  point  I 1  or  I2  to  the 
center  of  mass  of  the  grid  cell  in  which  this  point  lies,  compare  with  Figure  2.2b.  The 
value  d2  corresponds  to  the  distance  of  the  center  of  mass  of  an  h- box  to  the  center 
of  mass  of  a  grid  cell  in  situations  depicted  in  Figure  2.2a.  The  s  terms  describe 
limited  slopes  used  in  the  interpolation  process  as  described  below.  The  unlimited 
slope  corresponds  to  a  slope  obtained  by  linear  interpolation  between  neighboring  grid 
cells  in  the  x-direction,  i.e. 


si+|,j  —  (Qi+l,j  Qi,j)  /  ^x- 


Our  aim  is  to  find  appropriate  slope  limiters,  such  that  the  resulting  rotated  grid 
method  (with  upwind  fluxes)  avoids  oscillations  near  discontinuities. 

For  the  approximation  of  the  advection  equation  with  constant  advection  speed 
v  >  u  >  0  on  a  grid  with  Ax  =  Ay,  our  rotated  grid  finite  volume  method  takes  the 
form 


Q?r  =Q?J  -  ^Qlij  -  Qi- Id)  -  -  Q<j-i) 


=qZj  -  (Qij-i  +  Qij)  -  2  (Qi-u-i  +  Qi  hi) 

+  2  (5*+  Id-i  +  Si+bi)  dl  ~  2  (g*-|  J-i  +  Si- Id)  d\ 

At 


(2.17) 


Ax 


v[(l  -  d2)(Q™j  -  QZj^  +  d^QU,  ~  QZ-ij-i) 


Setting  all  the  s  terms  in  (2.17)  equal  to  zero  would  lead  to  the  monotone  method 
suggested  by  Roe  and  Sidilkover.  However,  our  main  goal  is  the  construction  of  a 
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second  order  rotated  grid  method  and  we  therefore  want  to  allow  slopes  that  are 
not  equal  to  zero  wherever  this  is  possible.  In  regions  where  the  solution  is  smooth 
neighboring  slopes  would  typically  have  about  the  same  size,  i.e.  we  have  in  particular 


2  ’•? 


'  2  ’•? 


and 


s*+| 


2  ->3 


'  2  ’•?  * 


(2.18) 


The  presence  of  such  slopes  in  (2.17)  would  not  destroy  the  monotonicity.  Near 
discontinuities  we  need  to  use  a  slope  limiter  to  enforce  a  relation  of  the  form  (2.18) 
to  hold  for  the  limited  slopes.  The  relation  (2.18)  motivates  the  use  of  a  two-sided 
slope  limiter.  This  can  be  seen  in  the  following  way:  In  grid  cell  (i,j)  the  relation 
ij  ps  Sj_i  j  suggests  using  a  left-sided  slope  limiter  in  order  to  limit  the  slope 
si+iy  This  can  for  instance  be  done  by  using  a  minmod  limiter  (see  for  instance 
[25])  that  can  be  written  in  the  form 


minmod(si+i)j,si_.|ij).  (2.19) 

In  grid  cell  (i  +  1  ,j)  the  limited  slope  should  satisfy  si+ij  ps  si+3j,  which  suggests 
using  a  right-sided  slope  limiter,  for  instance 


minmod(sj_|_i  ^ j).  (2.20) 

The  limited  slope  si+ij  is  used  to  calculate  the  upwind  point  needed  to  obtain  the 
numerical  flux  Fi+br  In  a  conservative  finite  volume  method  we  use  the  same  flux 
Fi+U  to  update  the  cell  average  of  the  conserved  quantity  in  grid  cell  ( i,j )  and 
( i  +  1,  j).  Therefore  we  need  a  combination  of  the  left-  and  right-sided  slope  limiter. 
This  can  be  obtained  by  applying  the  minmod  function  to  the  two  values  (2.19)  and 
(2.20)  which  results  in  the  two-sided  limiter 

si+iyj  -  minmod(si_i  j,  si+.y,si+3j).  (2.21) 

Our  numerical  experiments  show  that  it  is  indeed  necessary  to  consider  the  slopes 
from  both  sides  of  the  grid  cell  to  obtain  effective  slope  limiters.  The  next  test 
calculation  will  illustrate  this. 

Example  2.1.  We  consider  the  advection  equation  with  advection  speeds  u  —  0.8, 
v  —  1.  The  initial  values  contain  a  disk  with  radius  0.2  located  at  ( x,y )  =  (0.25,0.25) 
and  a  square  with  center  initially  located  at  (0.65,0.65). 

Figure  2.3b  shows  numerical  results  at  time  t  —  0.2  of  the  rotated  grid  method 
with  first  order  upwind  flux  without  the  use  of  slope  limiters  in  the  interpolation  step 
where  the  h- box  values  are  defined.  This  clearly  shows  oscillations  near  discontinuities. 
For  the  calculation  shown  in  (c)  we  used  a  one-sided  slope  limiter  of  the  form  (2.20). 
This  limiter  is  clearly  a  good  choice  for  the  test  problem  considered,  although  small 
oscillations  are  still  visible.  For  the  calculation  shown  in  (d)  we  used  a  one-sided 
limiter  of  the  form  (2.19),  which  turns  out  to  be  quite  ineffective.  In  (a)  we  show  the 
results  obtained  by  using  the  two-sided  limiter  (2.21).  The  same  limiting  process  can 
be  applied  in  more  general  situations,  e.g.  advective  transport  with  variable  velocity 
field  or  the  approximation  of  systems  of  conservation  laws. 

2.3.  Verification  of  the  high-resolution  rotated  grid  method  for  advec¬ 
tive  transport  problems.  In  order  to  approximate  the  2D  advection  equation  with 
a  rotated  grid  method,  we  choose  the  ^-direction  to  be  aligned  with  the  velocity  field. 
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Fig.  2.3.  Test  calculation  of  Example  2.1  with  different  slope  limiters.  All  calculations  use  a 
first  order  flux  calculation  but  linear  interpolation  to  determine  the  h-box  values.  400  x  400  grid 
cells,  CFLh  «  0.95,  t  =  0.2. 

The  fluxes  in  the  //-direction  are  then  equal  to  zero  and  we  do  not  have  to  construct 
h- boxes  in  this  orthogonal  direction.  At  each  cell  interface  we  approximate 

qt  +  \/u2  +  v'2qi  —  0, 

where  u  and  v  denote  the  advection  speed.  Let  F  be  approximated  by  the  Lax- 
Wendroff  flux  formula 

F  =  1-Vu2+v2(Ql  +  Qr )  -  ~{u2  +  v2){QR  -  Ql), 

where  QL  and  QR  are  the  /i-box  values  at  the  grid  cell  interface  considered  and  h  is 
the  length  of  the  h- box.  We  obtain  the  following  accuracy  result. 

Proposition  1.  Consider  the  approximation  of  sufficiently  smooth  solutions  of 
the  two-dimensional  advection  equation  by  a  rotated  grid  h-box  method.  The  direction 
of  the  h-boxes  are  aligned  with  the  velocity  field.  If  the  h-box  values  are  calculated  by 
(unlimited)  linear  interpolation  as  described  in  (2.16)  and  the  fluxes  F  are  approxi¬ 
mated  by  the  Lax-  Wendroff  method,  then  the  resulting  finite  volume  method  is  second 
order  accurate. 

We  omit  the  proof  which  is  based  on  Taylor  series  expansion  and  instead  show 
numerical  results  for  a  slightly  more  general  situation,  where  the  advection  speed  is 
locally  varying. 
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Example  2.2.  Solid  body  rotation.  The  time  independent  velocity  field  for 
solid  body  rotation  on  [0, 1]  x  [0, 1]  is  given  by 

u  —  —(y  —  0.5),  v  —  (x  —  0.5). 

As  initial  values  we  consider  a  smooth  hump  of  the  form 

q(x,  y,  0)  =  |  (1  +  cos  (tit (a;,  y))) , 

with 

r(x,  y)  =  min  ^ (x  -  x0)2  +  (y  -  y0)2,  r0)  /r0. 

The  parameters  are  set  to  Xo  —  0.25,  j/o  =  0.5,  r o  =  0.15.  At  time  t  —  2n,  i.e.  after 
one  rotation,  the  exact  solution  agrees  with  the  initial  values. 


Fig.  2.4.  Solution  of  the  rotating  flow  problem  after  one  rotation  using  high-resolution  flux 
formulas  on  a  grid  with  200  x  200  mesh  cells,  (a)  Contour  plot  of  solution,  (b)  comparison  of  exact 
and  numerical  solution  after  one  full  rotation  at  y  =  0.5. 


Figure  2.4  shows  numerical  results  for  the  rotating  flow  problem.  In  Figure  2.4a,  a 
contour  plot  of  the  numerical  solution  is  shown  after  one  full  rotation.  In  2.4b  we  show 
a  one-dimensional  slice  of  the  numerical  solution  at  y  —  0.5  together  with  the  exact 
solution  (solid  line).  Table  2.1  shows  the  Fi -errors  as  well  as  the  experimental  order  of 
convergence  (EOC)  after  one  rotation.  For  this  smooth  solution,  results  were  obtained 


(mx,  my)  /  EOC 

new  rotated  grid  method 
with  Lax-Wendroff  flux 

Roe  and  Sidilkover  method 
with  Lax-Wendroff  flux 

(100,100) 

0.27624 

0.33239 

(200,200) 

7.88745-2 

0.12841 

EOC 

1.81 

1.37 

(400,400) 

2.18381-2 

5.60114d-2 

EOC 

1.85 

1.20 

(800,800) 

5.83500d-3 

2.63392d-2 

EOC 

1.90 

1.09 

Table  2.1 

Accuracy  study  for  the  rotating  flow  problem. 
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using  our  new  rotated  grid  method  that  is  based  on  h- box  values  obtained  by  linear 
interpolation  as  described  above  and  a  Lax-Wendroff  flux  calculation.  The  numerical 
test  confirms  second  order  convergence.  We  also  show  the  error  and  numerical  order 
of  convergence  for  a  test  calculation  where  we  used  h- box  values  that  correspond  to 
the  optimal  positive  scheme  of  Roe  and  Sidilkover  combined  with  a  Lax-Wendroff  flux 
calculation.  As  expected  this  leads  to  a  first  order  accurate  approximation. 

2.4.  A  rotated  grid  method  for  systems  of  conservation  laws.  We  now 

extend  our  rotated  grid  method  to  the  approximation  of  two-dimensional  systems  of 
conservation  laws.  To  illustrate  the  performance  of  the  method  we  consider  the  Euler 
equations  of  gas  dynamics.  In  the  case  of  two  spatial  dimensions,  the  Euler  equations 
take  the  form 


pt  +  {pu)x  +  ( pv)y  —  0 
(. pu)t  +  {pu2  +p)x  +  ( puv)y  —  0 
{pv)t  +  ( puv)x  +  ( pv 2  +  p)y  -  0 
Et  +  (u(E  +  p))x  +  (V{E  +  p))y  —  0, 


where  (it,  v)  represents  the  velocity,  p  the  density,  p  the  pressure  and  E  the  energy. 
We  use  the  ideal  gas  equation  of  state 


E=~^r  +  \p{u2  +  v2). 
7  —  1  2 


The  approximation  of  two-dimensional  systems  of  conservation  laws  with  a  ro¬ 
tated  grid  method  requires  the  construction  of  h- boxes  in  the  £  and  17-direction.  The 
conserved  quantities  that  are  assigned  to  the  h- boxes  can  be  calculated  by  an  interpo¬ 
lation  procedure  analogous  to  that  of  the  advection  equation,  with  the  interpolation 
now  performed  for  a  vector  of  conserved  quantities.  Slope  limiters  for  each  conserved 
variable  are  used  to  prevent  unphysical  oscillations  near  discontinuities.  Once  the  h- 
box  values  are  known,  we  perform  the  flux  calculation  by  using  a  predictor- corrector 
method.  In  the  special  case  where  the  rotated  grid  is  aligned  with  the  Cartesian  grid 
coordinates,  the  predictor  corrector  method  agrees  with  Colella’s  multidimensional 
upwind  method  described  in  [10].  Using  a  predictor-corrector  method  the  numerical 
flux  F  (i.e.  the  numerical  flux  for  Equation  (2.4),  at  the  grid  cell  interface  (xi+i  ,Vj)) 
is  calculated  as 


=  F 


(<9f+ 1 , j (xi+i Mi tn+ 1 ) > Q?+ 1 ,3 ( xi+ § ,Vj ’tn+ 1 ) )  • 


(2.22) 


In  the  predictor  step  we  estimate  the  conserved  quantity  at  the  grid  cell  interface 
at  time  tn+ 2.  Those  values  are  then  used  in  the  corrector  step  (2.22).  Analogous 
predictor-corrector  steps  are  also  necessary  for  the  approximation  of  the  G  fluxes  at 
the  cell  interface  (xi+i,yj)  as  well  as  for  the  approximation  of  F  and  G  fluxes  at  cell 
interfaces  (xi,yj+i).  Using  Taylor  series  expansion,  we  find  that 


Qf+i  j(xi+h ,Vj,tn+i)  =  Q 


«+ih  _  nL’R 
i  + 1 J 


+ 

±-2"«e 


At 
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2 
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Here  A  is  the  Jacobian  matrix  of  the  flux  /  from  Equation  (2.8)  with  /  calculated 
by  (2.4).  The  distance  A£  is  equivalent  to  the  length  of  the  h- box  constructed  in 
^-direction.  The  predictor  step  can  be  performed  in  the  two  successive  steps 


Q 

Q 


L,R 

i+  Id 


=  Q 


L,R 

*+  Id 


L,R 

i+bi 


± 


A(_ma 


K 


At  d 


(2.23) 

(2.24) 


In  the  first  step  we  reconstruct  the  state  from  the  cell  center  to  the  cell  interface  by 
taking  into  account  only  the  characteristics  that  point  toward  the  interface.  This  is 
done  by  analogy  to  Colella’s  grid  aligned  method  described  in  [10].  In  the  second 
step  we  include  a  transverse  flux  difference.  This  step  is  necessary  in  order  to  obtain 
stability  for  time  steps  up  to  a  Courant  number  of  one.  (Note  that  in  a  rotated  grid 
method  for  the  advection  equation,  where  the  coordinate  direction  is  aligned  with 
the  characteristic  direction,  the  transverse  flux  component  would  be  equal  to  zero,  so 
we  did  not  have  to  include  a  transverse  flux  difference.)  In  order  to  incorporate  the 
transverse  flux  difference  into  the  predictor  step,  we  need  to  solve  Riemann  problems 
in  the  77-direction.  We  first  have  to  define  appropriate  initial  values  say  Q 1 ,  Q2  and  Q3. 
For  the  intermediate  value  we  use  our  h- box  value,  i.e.  Q2  —  Q'i  The  definition 

lT  2 

of  the  other  two  values  is  indicated  in  Figure  2.5.  From  the  center  of  mass  of  each 


Fig.  2.5. 
difference. 


Schematic  description  of  points  defined  for  the  calculation  of  the  transverse  flux 


h- box  we  find  points  along  a  line  in  the  transverse  direction  that  have  the  distance  h 
from  the  center  of  mass.  The  conserved  quantities  that  are  assigned  to  those  points 
are  obtained  by  interpolation  between  the  four  closest  neighboring  grid  cells  of  the 
underlying  Cartesian  grid.  Then  the  transverse  flux  difference  can  be  approximated 
by  (G(Q2,Q3)  ~  G{Q\ ,  C?  2 ) ) /h,  where  G  is  the  numerical  flux  obtained  by  solving 
Riemann  problems  for  Equation  (2.9). 

Below  we  perform  several  test  calculations,  using  model  problems  that  have  been 
studied  in  the  literature  with  grid  aligned  finite  volume  methods.  The  numerical 
results  obtained  with  our  rotated  grid  h- box  method  compare  well  with  those  obtained 
by  other  state  of  the  art  methods  for  hyperbolic  problems. 

Note  that  in  general  situations,  we  do  not  know  the  best  choice  for  the  direction 
of  the  rotated  grid  a  priori.  For  the  2D  Euler  equations,  different  approaches  may 
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be  possible.  A  simple  approach  is  to  choose  the  ^-direction  aligned  with  the  local 
velocity  field.  However,  if  we  want  to  approximate  shock  waves  that  are  not  aligned 
with  the  grid  a  better  choice  for  the  rotated  grid  might  be  to  choose  the  ^-direction 
aligned  with  the  pressure  gradient.  With  such  an  approach  we  hope  to  get  sharper 
approximation  of  shock  waves.  By  comparing  first  order  accurate  methods  a  much 
sharper  approximation  of  shock  waves  that  are  not  aligned  with  the  grid  can  indeed 
be  obtained  by  using  a  rotated  grid  method  that  is  aligned  with  the  shock  instead  of 
a  grid  aligned  method,  see  [27].  However,  high-resolution  grid  aligned  methods  lead 
to  very  accurate  approximation  of  shock  waves  that  are  not  aligned  with  the  grid  and 
we  do  not  see  any  further  improvement  by  using  our  rotated  grid  method.  Instead  we 
will  see  that  our  rotated  grid  method  compares  well  with  grid  aligned  high-resolution 
methods. 

We  first  demonstrate  the  performance  of  our  high-resolution  rotated  grid  method 
for  a  two-dimensional  Riemann  problem,  as  suggested  in  [35].  Note  that  this  test  case 
was  also  considered  in  [28],  where  several  numerical  methods  were  compared. 

Example  2.3.  We  approximate  the  solution  of  the  two-dimensional  Euler  equa¬ 
tions  with  piecewise  constant  initial  values.  Let  the  index  indicate  the  quadrant,  then 
the  initial  values  have  the  following  form,  pi  —  0.4,  p\  —  0.5313,  rti  =  0.0,  Vi  — 
0.0,  p2  —  1.0,  p2  —  1.0,  u2  —  0.7276,  v2  —  0.0,  p3  —  1.0,  pz  —  0.8,  u3  —  0.0,  v3  — 
0.0, pi  —  1.0, p4  —  1.0,  «4  =  0.0, 14  =  0.7276.  The  ratio  of  specific  heat  is  set  to 
7  =  1.4  and  the  computational  domain  is  [0,1]  x  [0,1].  All  neighboring  Riemann 
problems  resulting  from  this  initial  value  problem  lead  to  shock  wave  solutions. 

Figure  2.6  shows  results  of  a  numerical  approximation  of  density  and  pressure.  To 


Density  at  time  t=0.25  Pressure  at  time  t=0.25 


Fig.  2.6.  Contour  plots  of  density  and  pressure  for  Example  2.3.  400  x  400  grid  cells.  The 
direction  of  the  rotated  grid  is  determined  from  the  pressure  gradient. 


make  comparison  with  the  schemes  shown  in  [28]  possible,  we  use  the  same  resolution, 
i.e.  a  grid  with  400  x  400  mesh  cells.  With  our  rotated  grid  method,  near  shock 
waves  the  grid  direction  is  aligned  with  the  pressure  gradient  while  the  direction  of 
the  underlying  Cartesian  grid  is  used  away  from  discontinuities.  We  get  a  sharp 
approximation  of  the  shock  structure  and  also  demonstrate  that  our  limiting  avoids 
unphysical  oscillations. 

Next  we  perform  an  accuracy  study  by  comparing  the  results  of  a  2D  radially 
symmetric  problem  with  the  numerical  solution  of  a  one  dimensional  calculation  with 
geometric  source  term  that  models  the  radial  symmetry. 
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Example  2.4.  We  consider  the  2D  Euler  equations  with  initial  values  that  have 
zero  velocity.  Initial  values  for  density  and  energy  are  given  by  the  smooth  function 


P{x,y,  0)  =  E(x,y, 


1  —  0.1(cos(47rr)  —  1) 
1 


0  <  r  <  0.5 
r  >  0.5, 


with  r  —  \J x2  +  y2.  The  computational  domain  is  [0, 1]  x  [0, 1],  with  solid  wall  bound¬ 
ary  condition  at  x  —  0  and  y  —  0. 

Numerical  results  for  this  test  problem  are  shown  in  Figure  2.7.  An  experimental 


Density  at  time  t=0.5  Scatter  plot  of  density  at  time  t=0.5 


Fig.  2.7.  Solutions  of  Example  2.4  on  a  grid  with  100  x  100  mesh  cells,  (a)  Contour  plot  of 
density,  (b)  Scatter  plot  of  2D  solution  together  with  ID  reference  solution  (solid  line).  Second  order 
method  without  limiters.  The  £- direction  of  the  rotated  grid  is  aligned  with  the  radial  direction. 


convergence  study  is  shown  in  Table  2.2,  which  confirms  that  our  method  is  second 
order  accurate  also  for  this  nonlinear  system.  Now  the  radial  direction  describes  the 
^-direction  for  our  rotated  method. 


(mx,  my)  /  EOC 

P 

pu 

E 

(50,50) 

1.9630d-4 

1.3794d-4 

2.7403d-4 

(100,100) 

5.0255d-5 

3.5021d-5 

7.0454d-5 

EOC 

1.96 

1.98 

1.96 

(200,200) 

1.3376d-5 

9.2185d-6 

1.8726d-5 

EOC 

1.91 

1.93 

1.91 

Table  2.2 


Accuracy  study  for  Example  2.f  using  the  second  order  rotated  grid  method  without  limiters.  A 
rotated  grid  method  was  used  in  which  the  £- direction  is  aligned  with  the  radial  direction.  Errors  are 
calculated  in  the  L\  norm.  The  ID  reference  solution  was  calculated  on  a  grid  with  40000  mesh  cells. 
For  the  2D  simulation  we  used  a  discretization  with  A t  =  Ax  =  A y,  which  leads  to  CFL  R*  0.8. 


This  test  case  was  also  considered  in  [24],  where  results  for  the  high-resolution 
wave  propagation  algorithm  were  shown.  Our  results  are  approximately  equivalent 
to  the  results  documented  in  [24],  The  wave  propagation  method  or  Colella’s  grid 
aligned  method  [10]  as  well  as  other  grid  aligned  methods  are  of  course  less  expensive 
than  our  rotated  grid  method.  We  need  to  solve  about  twice  as  many  Riemann 
problems  as  these  other  methods.  However,  our  main  motivation  for  the  development 
of  the  high-resolution  rotated  grid  method  is  to  overcome  the  small  cell  problem  for 
Cartesian  grid  methods  with  embedded  boundaries.  In  that  case,  we  only  want  to  use 
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the  rotated  grid  method  near  the  embedded  boundary,  while  a  classical  grid  aligned 
method  can  be  used  away  from  the  boundary. 

3.  The  embedded  boundary  method.  Based  on  our  high-resolution  rotated 
grid  method  we  now  describe  a  high-resolution  rotated  grid  embedded  boundary 
method.  The  resulting  method  is  a  refinement  of  the  embedded  boundary  method  of 
Berger  and  LeVeque  [4],  [5].  We  assume  that  the  embedded  boundary  is  described 
by  a  piecewise  linear  segment  that  cuts  through  each  boundary  cell  once.  Then  each 
irregular  grid  cell  is  a  polygon  with  at  most  five  sides.  Thus  a  finite  volume  method 
has  the  form 


Qu1  =  Qli- 


A  t 


ei+iFi+ij 


L  F, 


h,i  +  ei+iGi,j+h 


Here  F  and  G  denote  the  fluxes  in  the  x  and  ^-direction  respectively.  II  is  the  flux 
per  unit  time  through  the  irregular  side  of  the  cell,  the  solid  wall  boundary.  The  value 
atj  describes  the  size  of  the  grid  cell  and  £i+i  is  the  length  of  the  interface  between 
cell  (i,j)  and  (i  +  l,j).  Note  that  £  may  be  zero  if  the  corresponding  grid  cell  interface 
is  completely  covered  by  the  boundary.  In  order  to  obtain  a  conservative  boundary 
flux,  we  construct  h- boxes  normal  to  the  boundary.  This  is  illustrated  in  Figure  3.1a. 
To  allow  time  steps  based  on  a  CFL  condition  that  is  appropriate  for  the  regular  part 
of  the  grid,  the  domain  of  dependence  has  to  be  large  enough.  By  analogy  with  our 
rotated  grid  method  for  regular  Cartesian  grids,  we  use  h- boxes  of  length  described 
in  Equation  (2.11). 

At  the  embedded  boundary  the  h- box  method  retains  stability  by  constructing  a 
finite  volume  scheme  where  the  flux  differences  in  each  cut  cell  are  of  the  order  of  the 
size  of  the  small  cell.  This  balances  the  term  a,  j  in  the  denominator.  The  direction 
normal  to  the  boundary  segment  is  the  77-direction  of  our  rotated  grid  method.  For  a 
small  triangular  shaped  cell  the  h- boxes  constructed  in  the  77-direction  are  shown  in 
Figure  3.1a,  b,  c.  Those  h- boxes  overlap  except  on  an  area  of  the  size  of  the  small  grid 
cell.  In  Figure  3. Id  and  e,  we  show  the  h- boxes  constructed  in  orthogonal  direction 
(the  ^-direction  of  our  rotated  grid  method).  Those  h- boxes  also  overlap  except  on 
an  area  of  the  size  of  the  small  cell.  Since  the  flux  calculation  is  Lipschitz  continuous, 
this  gives  us  the  cancellation  property  needed  to  obtain  stability  in  the  small  cell.  A 
construction  of  the  boundary  h- boxes  in  direction  normal  to  the  boundary  allows  us 
to  construct  a  conservative  method. 

In  the  two-dimensional  case  we  do  not  have  a  stability  proof  of  our  method  at  cut 
cells.  However,  the  heuristic  argument  given  above  is  a  generalization  of  the  stability 
argument  developed  in  the  one-dimensional  case.  We  have  proved  stability  of  the 
second  order  accurate  h- box  method  on  one-dimensional  irregular  grids  in  [2]. 

In  order  to  obtain  high-accuracy  near  the  boundary  of  an  embedded  object,  we 
have  to  find  accurate  interpolation  formulas  for  the  calculation  of  the  h- box  values. 
Furthermore,  these  h- box  values  have  to  be  constructed  in  a  way  that  avoids  oscil¬ 
lations  near  discontinuities.  This  can  be  obtained  by  a  generalization  of  the  limiting 
described  in  section  2.2. 

In  order  to  measure  the  accuracy  of  our  embedded  boundary  method  at  the 
boundary  as  well  as  inside  the  domain,  we  use  two  different  norms.  Inside  the  domain 
we  can  calculate  the  error  with  a  standard  L 1  -norm  of  the  form 


E ij  \Qi,j  -  <i(xi,yj)\ai,j 

Eij  \q(xuyj)\aij 
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Fig.  3.1.  (a)  Construction  of  a  boundary  h-box,  i.e.  in  rj-direction  of  our  rotated  grid,  (b),  (c) 
Construction  of  h-boxes  in  rj-direction.  (d),  (e)  Construction  of  h- boxes  in  direction. 


Here  Qij  is  the  numerical  solution  in  grid  cell  (i,j)  and  q(xt,  i/j)  is  the  exact  solution 
at  the  center  of  mass  of  this  grid  cell.  The  area  of  the  grid  cell  (i,j)  is  a,_j.  To  measure 
the  boundary  error  we  sum  over  all  cut  cells  and  calculate  the  norm 

_  T,(i,j)£K\Qi,i  -  q(xhyj)\hj 
b~  Z^eK 

where  bij  is  the  length  of  the  boundary  segment  in  grid  cell  K  is  the  set  of  all 

cut  cells  along  the  embedded  boundary. 

A  second  order  accurate  approximation  of  the  solution  over  the  whole  domain 
can  be  obtained,  even  if  the  boundary  is  only  approximated  with  first  order  accuracy. 
This  is  due  to  the  lower  dimensionality  of  the  boundary  segment,  see  Gustafsson  [17]. 
However,  in  several  applications  one  may  be  especially  interested  in  the  solution  near 
the  boundary.  In  such  cases  an  accurate  boundary  treatment  is  essential.  We  start 
the  description  of  our  embedded  boundary  method  by  presenting  an  approach  for 
the  advection  equation  with  variable  velocity  field.  The  velocity  field  is  described 
by  a  stream  function  that  may  be  given  by  an  analytical  formula  or  by  a  numerical 
approximation.  By  studying  advective  transport  problems  first,  we  again  only  have 
to  construct  h- boxes  for  our  rotated  grid  method  in  the  direction  aligned  with  the 
velocity  field. 

3.1.  Advection  equation.  In  this  subsection  we  construct  an  embedded  bound¬ 
ary  rotated  grid  method  for  the  approximation  of  advective  transport  problems  with 
variable  velocity  field.  We  assume  that  the  velocity  field  near  the  boundary  is  tangen¬ 
tial  to  the  boundary  and  may  be  given  by  a  stream  function.  Such  a  situation  may 
for  instance  describe  the  advection  of  a  tracer  through  a  field  with  irregular  inclusions 
to  model  porous  media  flow.  This  situation  can  be  simulated  by  a  rotated  grid  h- box 
method  with  h- boxes  in  the  ^-direction  only. 

By  analogy  with  our  rotated  grid  h- box  method  for  regular  Cartesian  grids,  we 
want  to  construct  h- box  values  that  represent  a  linear  interpolation  of  the  conserved 
quantities  to  the  center  of  mass  of  the  h- boxes.  We  first  calculate  for  each  h- box  the 
area  fractions  of  all  polygons  that  represent  the  intersection  of  the  h- box  with  the 
underlying  Cartesian  grid  cells  that  are  overlapped  by  the  h- box.  As  is  clear  from 
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Figure  3. Id  and  e,  an  h- box  in  the  ^-direction  intersects  at  most  two  Cartesian  grid 
cells.  We  calculate  the  center  of  mass  of  these  two  polygons.  We  also  need  to  know 
the  center  of  mass  of  the  underlying  grid  cells  that  are  overlapped  by  the  h- box  as 
well  as  a  gradient  of  the  conserved  quantities  in  each  grid  cell.  Using  this  information 
we  can  interpolate  the  conserved  quantities  at  the  center  of  mass  of  the  polygons  that 
represent  the  intersection  of  the  h- box  with  the  grid  cells.  The  conserved  quantity 
that  represents  an  interpolation  at  the  midpoint  of  the  h- box  can  then  be  obtained  as 
an  area- weighted  sum  of  the  interpolated  values  at  the  center  of  mass  of  the  polygon. 
In  order  to  calculate  the  gradients  of  the  conserved  quantities  in  grid  cells  near  the 
boundary,  we  use  a  least  square  procedure.  A  limiter  is  used  to  prevent  oscillations 
near  discontinuities. 

To  test  our  embedded  boundary  method  we  first  consider  advective  transport 
inside  an  annulus. 

Example  3.1.  We  consider  solid  body  rotation  inside  an  annulus  defined  by  two 
concentric  circles.  The  radius  of  the  inner  circle  is  Ri  —  0.75  and  the  radius  of  the 
outer  circle  is  i?2  =  1.25.  The  annulus  is  embedded  into  a  Cartesian  grid  of  the  size 
[—1.5, 1.5]  x  [—1.5, 1.5].  We  use  initial  values  of  the  form 

q{x,y,t)  =  w(6  -  tt/2) 


with 


w{6)  —  0.5 


erf 


7t/6  —  6  ^ 

7W o) 


+  erf 


7t/6  +  6  ^  ^ 

Tt/foOylj  ' 


The  velocity  field  describing  solid  body  rotation  was  obtained  using  the  stream  function 

ip(x,y)  =  0.27r(i?2  -  r2), 

with  r  —  (x2  + 1/2)1/2.  One  rotation  is  completed  at  time  t  —  5. 

Figure  3.2  shows  contour  plots  of  the  numerical  solution  at  different  time  steps. 
The  solution  at  the  two  boundaries  is  shown  in  Figure  3.3.  We  used  our  second  order 


Initial  values 


Solution  after  half  rotation  Solution  after  one  rotation 


Fig.  3.2.  Solution  of  the  annulus  test  problem  on  a  grid  with  800  x  800  grid  cells. 

rotated  grid  method  for  regular  grid  cells  and  the  rotated  grid  embedded  boundary 
method  along  the  boundary.  For  the  fine  discretization  used  here,  with  800  x  800  grid 
cells  on  the  rectangular  computational  domain,  we  have  many  very  small  cut  cells, 
the  smallest  of  which  is  about  10  5  times  the  size  of  a  regular  grid  cell.  However, 
this  does  not  cause  any  stability  or  accuracy  problem.  In  Table  3.1  we  show  results 
of  a  convergence  study,  comparing  the  solution  inside  the  domain  as  well  as  along 
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Solution  along  inner  boundary  Solution  along  outer  boundary 


Fig.  3.3.  Solution  of  the  annulus  test  problem  along  the  inner  and  outer  boundary  after  one 
rotation  for  0  <  6  <  tv. 


the  boundaries  with  the  exact  solution.  One  can  see  that  the  highest  accuracy  was 
obtained  along  the  outer  boundary,  while  the  lowest  accuracy  was  obtained  along  the 
inner  boundary.  This  may  not  be  surprising,  since  the  solution  profile  is  resolved  with 
more  points  per  transition  layer  as  the  radius  increases. 


mx  x  my  /  EOC 

domain 

outer  boundary 

inner  boundary 

100  x  100 

7.50839d-2 

3.87471d-2 

0.13909 

200  x  200 

2.03198d-2 

1.01406d-2 

4.01216d-2 

EOC 

1.88 

1.93 

1.79 

400  x  400 

5.14101d-3 

2.56584d-3 

1.02206d-2 

EOC 

1.98 

1.98 

1.97 

800  x  800 

1.29810d-3 

6.45543d-4 

2.56709d-3 

EOC 

1.99 

1.99 

1.99 

Table  3.1 

Convergence  study  for  annulus  test  problem. 


Note  that  advective  transport  inside  an  annulus  was  also  studied  in  Calhoun  and 
LeVeque  [7],  where  a  Cartesian  grid  embedded  boundary  method  for  the  advection 
diffusion  equation  was  developed.  Their  method  remains  stable  in  the  case  of  pure 
advection.  However,  in  the  advection  case  the  accuracy  at  the  boundary  was  reduced, 
which  caused  an  unphysical  boundary  layer. 

Numerical  convergence  studies  can  also  be  found  in  Pember  et  al.  [32]  as  well  in 
Forrer  and  Jeltsch  [16].  Those  tests  document  a  loss  of  accuracy  at  the  boundary  of 
those  earlier  Cartesian  grid  embedded  boundary  methods. 

3.2.  Euler  equations.  In  order  to  obtain  an  embedded  boundary  method  for 
the  Euler  equations  we  construct  in  addition  to  the  h- boxes  in  the  ^-direction  also 
h- boxes  in  the  77-direction.  At  the  boundary  segment,  which  is  a  straight  line  for  each 
cut  cell,  we  construct  a  boundary  h- box  normal  to  the  boundary.  The  length  of  this 
h- box  is  again  determined  by  Equation  (2.11),  where  (it, v)  may  now  represent  the 
normal  vector  at  the  boundary.  As  depicted  in  Figure  3.1a,  we  construct  a  boundary 
h- box  that  points  into  the  flow  domain  (the  inbox)  as  well  as  one  that  points  into 
the  boundary  (the  outbox).  The  conserved  quantities  assigned  to  the  outbox  are  the 
reflected  values  of  the  conserved  quantities  of  the  inbox.  For  the  Euler  equations  this 
means  that  for  a  vector  of  conserved  quantities  (given  in  Cartesian  coordinates)  we 
first  determine  (by  rotation)  the  normal  and  tangential  velocity  in  (£,  77)-coordinates. 
We  denote  those  velocities  by  u  and  v.  The  outbox  values  are  then  obtained  from  the 


20 


C.  HELZEL,  M.J.  BERGER  AND  R.J.  LEVEQUE 


inbox  values  by  setting 


P outbox  —  Pinboxt  ^outbox 


^inbox  i  ^ outbox  —  ^inboxiP outbox  — Pinbox 


This  is  a  widely  used  procedure  to  obtain  boundary  fluxes  that  simulate  a  reflecting 
boundary,  see  for  instance  [25].  It  allows  us  to  calculate  the  boundary  flux  by  using 
the  same  Riemann  solver  routines  that  are  used  inside  the  computational  domain.  In 
a  typical  situation,  where  a  reflecting  boundary  is  aligned  with  a  grid  cell  interface, 
the  outbox  would  be  called  a  ghost  cell.  The  conserved  quantities  assigned  to  the 
inbox  are  obtained  in  the  same  way  as  described  before  for  the  h- boxes  in  the  In¬ 
direction.  For  each  grid  cell  that  is  overlapped  by  the  inbox  we  calculate  a  polygon 
that  represents  the  intersection  of  the  inbox  with  this  grid  cell.  We  then  reconstruct 
the  conserved  quantities  at  the  center  of  mass  of  each  of  these  polygons  and  calculate 
an  area-weighted  average,  which  is  equivalent  to  a  piecewise  linear  interpolation  of 
the  conserved  quantities  at  the  center  of  mass  of  the  inbox. 

It  remains  to  calculate  the  h- box  values  in  the  77-direction  at  the  cell  interfaces 
near  the  boundary,  as  indicated  in  Figure  3.1b  and  c.  The  calculation  of  the  conserved 
quantities  assigned  to  the  h- boxes  that  point  into  the  flow  domain  follows  the  proce¬ 
dure  used  for  the  calculation  of  the  values  assigned  to  the  inbox  as  well  as  the  h- boxes 
in  the  ^-direction.  The  conserved  quantities  assigned  to  the  h- box  that  points  toward 
the  boundary  is  slightly  more  complicated.  A  fraction  of  this  h- box  can  be  inside  the 
boundary.  We  first  calculate  the  part  of  the  h- box  that  intersects  with  the  boundary 
and  lies  outside  the  flow  domain.  This  polygon  is  then  reflected  at  the  boundary 
segment.  For  the  reflected  polygon  we  can  now  calculate  the  intersection  with  grid 
cells  that  are  overlapped  by  this  polygon.  We  then  calculate  a  linear  interpolation  of 
the  conserved  quantities  at  the  center  of  mass  of  the  reflected  polygon.  This  vector  of 
conserved  quantities  is  reflected  and  assigned  to  the  part  of  the  h- box  that  lies  inside 
the  boundary.  For  the  fraction  of  the  h- box  that  lies  inside  the  flow  domain  we  also 
calculate  a  linear  interpolation  of  the  conserved  quantities  at  the  center  mass  of  this 
polygon.  The  area  weighted  average  of  conserved  quantities  that  are  assigned  to  those 
two  polygons  is  used  as  the  h- box  value.  Note  that  depending  on  the  geometry  of 
the  embedded  boundary  a  part  of  an  h- box  constructed  in  the  ^-direction  might  also 
intersect  with  the  boundary,  in  this  case  we  proceed  as  described  above  for  h- boxes 
in  the  77-direction. 

Once  the  h- box  values  are  defined,  we  calculate  fluxes  at  the  cell  interfaces  of  a 
cut  cell  using  the  same  rotated  grid  predictor-corrector  method  as  outlined  in  section 
2.4.  Note  that  this  also  includes  the  calculation  of  transverse  flux  differences.  In  order 
to  calculate  the  transverse  flux  differences  we  define  vectors  Q 1,  Q2  and  Q3  of  the 
conserved  quantities,  compare  with  section  2.4.  Here  Q 2  is  again  the  h- box  value. 
The  other  two  points  lie  in  the  orthogonal  direction  a  distance  h  away  from  the  center 
of  mass  of  the  h-hox.  If  those  two  points  are  inside  the  computational  domain,  then 
linear  interpolation  between  neighboring  grid  cells  is  used  to  assign  the  values  Q 1  and 
Q 3.  If  one  or  both  of  the  points  that  are  constructed  in  the  orthogonal  direction  are 
outside  the  domain,  then  we  first  calculate  the  reflected  point  (which  is  inside  the  flow 
domain).  We  assign  a  vector  of  conserved  quantities  to  the  reflected  point.  The  vector 
used  in  the  calculation  of  the  transverse  flux  difference  can  be  obtained  by  reflection 
of  the  conserved  quantities. 

Below  the  performance  of  our  embedded  boundary  method  is  demonstrated  for 
shock  reflection  from  a  wedge  as  well  as  shock  reflection  from  a  cylinder. 


A  ROTATED  GRID  METHOD  FOR  EMBEDDED  GEOMETRIES 


21 


Shock  reflection  from  a  wedge.  Here  we  study  the  approximation  of  a  Mach 
10  shock  wave  reflected  from  a  30  degree  wedge.  The  state  in  front  of  the  shock  is 


Entropy  at  time  t=0.2 


JL 

0.5  1  1.5  2  2.5  3 

Fig.  3.4.  Shock  reflection  from  a  30  degree  wedge.  Parameter  values:  mx=  450,  my  =  300, 
size  of  the  smallest  cut  cell  4.2  x  10-6;  high-resolution  rotated  grid  method  with  minmod  limiter, 
CFL  <  1 .  Away  from  the  boundary  the  coordinate  direction  is  aligned  with  the  grid  and  no  rotation 
is  used. 

p  —  1.4,  u  —  v  —  0,  p  —  1.  At  the  embedded  boundary  we  use  a  rotated  grid  which 
is  aligned  with  the  boundary.  Inside  the  domain  we  use  the  grid-aligned  multidimen¬ 
sional  upwind  method  of  Colella.  For  this  test  case  the  size  of  the  smallest  cut  cell 
was  4.2  x  10  6  times  the  size  of  a  regular  grid  cell.  This  did  not  cause  any  stability 
problem  and  moreover  we  get  an  accurate  approximation  near  the  boundary  as  well 
as  in  the  whole  domain.  Figure  3.4  shows  numerical  results  for  density  and  entropy. 

Shock  reflection  from  a  cylinder.  Finally  we  consider  the  performance  of 
our  embedded  boundary  method  for  the  Euler  equations  in  the  presence  of  a  curved 
boundary.  As  a  numerical  test  calculation  we  consider  shock  reflection  from  a  circular 
cylinder.  Initially  a  Mach  2  shock  is  located  at  x  —  —0.3.  The  state  in  front  of  the 
shock  is  again  p  —  1.4,  u  —  v  —  0,  p  —  1.  Figure  3.5a,  b  shows  contour  plots  for 
pressure  and  density  obtained  by  our  Cartesian  grid  embedded  boundary  method. 

We  determine  the  direction  of  the  rotated  grid  at  each  cut  cell  interface  such  that 
the  ^-direction  is  tangential  to  the  boundary.  A  smooth  transition  from  the  rotated 
grid  to  the  mesh-aligned  grid  for  grid  cells  near  the  boundary  was  used.  In  Figure 
3.5c,  d,  we  show  results  for  the  same  test  problem  obtained  by  the  curvilinear  version 
of  CLAWPACK  [19],  see  LeVeque  [25]  for  a  description  of  the  algorithm.  This  uses  a 
logically  rectangular  grid  with  90  x  360  grid  cells,  i.e.  the  inner  boundary  is  resolved 
by  360  full  cells.  Although  a  curvilinear  grid  is  of  course  a  very  good  choice  for 
this  simple  situation,  we  can  see  from  this  comparison  that  our  embedded  boundary 
method  also  leads  to  an  accurate  approximation  of  the  solution  structure. 

In  order  to  further  study  the  accuracy  of  the  embedded  boundary  method,  we 
also  compare  the  solution  along  the  boundary,  i.e.  the  solution  in  the  cut  cells  ob¬ 
tained  by  our  embedded  boundary  method  and  the  solution  along  the  first  row  of 
grid  cells  for  the  curvilinear  grid  cells.  This  comparison  is  shown  in  Figure  3.6  for 
a  quite  coarse  grid,  i.e.  60  cells  along  the  boundary  (cut  or  full  respectively).  Note 
that  in  order  to  make  comparisons  easier,  the  grids  are  much  coarser  than  the  ones 
used  for  the  calculation  shown  in  Figure  3.5.  Again  we  can  see  that  our  embedded 
boundary  method  leads  to  an  accurate  approximation.  It  is  noticeable  however,  that 
the  curvilinear  grid  method  leads  to  a  sharper  approximation  of  the  shock  wave.  Such 
a  behavior  was  already  observed  for  our  one-dimensional  irregular  grid  h- box  method 


Density  at  time  t=0.2 
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Density  at  time  t=0.25  Pressure  at  time  t=0.25 


Fig.  3.5.  Reflection  of  a  Mach  2  shock  wave  from  a  circular  cylinder,  (a),  (b)  Density  and 
pressure  obtained  on  Cartesian  grid  with  300  x  300  grid  cells,  364  cu t  ce/Zs.  (c),  (d)  Density  and 
pressure  for  curvilinear  grid  method  on  a  90  x  360  mesh.  Note  that  on  a  curvilinear  grid  the  shock 
is  not  aligned  with  the  mesh.  Therefore,  simple  extrapolating  boundary  conditions  (used  here)  can 
cause  some  weak  waves  that  are  moving  back  into  the  computational  domain. 


presented  in  [2].  Since  several  h- boxes  may  overlap  a  single  discontinuity,  the  method 
can  have  more  smearing  than  usually  observed  on  a  regular  Cartesian  grid. 


Density  along  the  boundary  —  Cartesian  grid  method  Density  along  the  boundary  —  Curvilinear  grid  method 


Fig.  3.6.  Density  along  the  boundary  for  reflection  of  Mach  2  shock  from  a  cylinder.  Solution 
is  shown  for  time  t  =  0.2.  (a)  Cartesian  grid  method  on  a  50  x  50  grid  with  60  cut  cells  along 
the  boundary;  (b)  curvilinear  grid  with  60  regular  grid  cells  along  the  boundary.  The  solid  line  is  a 
reference  solution  obtained  by  a  curvilinear  grid  calculation  on  a  fine  mesh. 
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As  soon  as  the  geometry  becomes  more  complex,  the  construction  of  a  body  fitted 
logically  rectangular  mesh  becomes  more  difficult.  This  is  for  instance  the  case  in  our 
last  test  problem,  where  we  consider  Mach  reflection  from  two  circular  cylinders.  The 
Cartesian  grid  embedded  boundary  method  can  of  course  handle  such  a  situation. 
Contour  plots  of  the  pressure  obtained  with  our  Cartesian  grid  embedded  boundary 
method  are  shown  in  Figure  3.7. 


Pressure  -  Cartesian  grid  method  Pressure  -  Cartesian  grid  method 


Fig.  3.7.  Reflection  of  a  Mach  3  shock  wave  from  two  circular  cylinders.  Contour  plots  of 
pressure  at  two  different  time  steps  (t=0.096,  t=0.16)  calculated  on  an  underlying  Cartesian  grid 
with  300  x  300  mesh  cells.  Initially  a  Mach  3  shock  was  located  at  x  =  0.2. 


There  has  been  recent  work  on  the  general  problem  of  implementing  boundary 
conditions  at  curved  boundaries,  for  example  see  Dadone  and  Grossman  [11],  [12]  for 
modifications  of  the  usual  technique  for  using  ghost  cells  for  reflecting  boundaries 
cells  for  steady  state  problems.  We  expect  that  similar  improvements  for  boundary 
conditions  at  curved  reflecting  boundaries  could  be  incorporated  into  our  h- box  imple¬ 
mentation  at  curved  boundaries  as  well.  We  plan  to  investigate  this  further  for  both 
curvilinear  grids  as  well  as  for  our  Cartesian  grid  embedded  boundary  method.  A 
higher  order  representation  of  the  geometry  may  also  be  necessary  in  order  to  further 
increase  the  accuracy,  see  also  [1]. 

4.  Conclusions.  We  developed  a  high-resolution  rotated  grid  h- box  method, 
which  is  second  order  accurate  for  smooth  solutions  of  conservation  laws  and  pre¬ 
vents  unphysical  oscillations  near  discontinuities.  This  method  was  incorporated  in  a 
Cartesian  grid  embedded  boundary  method.  The  resulting  method  leads  to  a  stable, 
accurate  and  conservative  approximation  of  conservation  laws  in  the  presence  of  ir¬ 
regular  embedded  boundaries.  Our  new  rotated  grid  h- box  method  is  more  accurate 
than  other  existing  Cartesian  grid  embedded  boundary  methods  for  conservation  laws. 
In  order  to  further  increase  the  accuracy  at  the  boundary,  it  would  be  interesting  to 
study  boundary  treatments  that  use  a  higher-order  representation  of  the  boundary. 
In  our  current  approach  the  boundary  segment  is  approximated  piecewise  linearly  in 
each  cut  cell.  Higher  accuracy  with  the  current  approach  can  also  be  obtained  by 
adaptive  mesh  refinement.  By  using  different  patches  of  Cartesian  grids  (as  for  in¬ 
stance  described  in  [3]),  we  can  apply  our  embedded  boundary  treatment  on  every 
level  of  refinement  in  exactly  the  same  form  as  described  in  this  paper  for  a  single 
grid. 
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